% Computes impulse response of states and jumps from Klein solution.
%    then calls computeirf.m to calculate IRFs of all other variables.
%
% Jumps are: V, C, p. 
% States are: z, Phihat.  
%      z is an exogenous state with shocks; Phihat is a predetermined, endogenous state.
% All state and jump variables represent deviations from the steady state,
%      so STATEDYNAMICS must have all eigenvalues <= 1.
%
% Klein solution is:
%          states_t+1 = STATEDYNAMICS*states_t + [shocks_t+1 ; 0]
%             jumps_t = JUMPS*states_t.
%
% for simplicity, 'time' is same as index of row vectors: starts at time=1.

% SET TIME 1 DISTRIBUTION (both dist and shock can be 0 or nonzero at initial time)
fprintf('\nDISTSIM: STARTING FROM STEADY STATE DISTRIBUTION')
fprintf('\n\n') 
shiftdist = Pdist;

initdist = shiftdist - Pdist;   % deviation from steady state
initdist = initdist(:);

initscalarstates = zeros(nss,1); % scalar states start from steady-state 

% PUT TOGETHER ALL TIME 1 INITIAL CONDITIONS
initSTATE = [time1Rshock; time1TFPshock; initdist; initscalarstates];   % TIME 1 STATE
initJUMPS = JUMPS*initSTATE;

STATEHISTORY = NaN*zeros(nz+nPhi+nss,TT);
JumpHistory = NaN*zeros(nV+nsj,TT);

STATEHISTORY(:,1) = initSTATE;           % time 1 state
JumpHistory(:,1) = initJUMPS;            % time 1 jumps

% initialize time 1 state for loop:
STATENOW = initSTATE;

% SPECIFY ZERO SHOCKS for the distribution
zeroshocks = zeros(nump+nss,1);  

% NOW SIMULATE DYNAMICS
for time = 2:TT
   STATENOW = STATEDYNAMICS*STATENOW + [Rshocks(time); TFPshocks(time); zeroshocks];
   JumpHistory(:,time)  = JUMPS*STATENOW;
   STATEHISTORY(:,time) = STATENOW;
end




